clear all
set more off
set scheme s2color

// Historical balancing tests

global input "...\replication\data"
global output "...\replication\outputs"

use "$input\com_additional_analysis.dta"

****************** Get other vars ******************
merge 1:1 insee using "$input\com_balance.dta"
drop if _merge != 3
drop _merge
merge 1:1 insee using "$input\france_cassini_cities.dta"
drop if _merge != 3
drop _merge

gen inter_roman = (OBJECTID_1 != 0)
gen inter_diocese = (FID_dioces != 0)
gen cassini_road = (CLASS != "")
gen cassini_city = (NAME != "")

ren insee com

****************** Costmetics ******************
** Remove Corsica
drop if substr(com,2,1) == "A"
drop if substr(com,2,1) == "B"
destring com code_dept, replace


****************** Get dist border ******************
merge 1:1 com using "$input\com_arcgis.dta"
drop if _merge != 3
drop _merge
****************** Get dist instrumental border ******************
merge m:1 com using "$input\com_dist_instr.dta"
drop if _merge != 3
drop _merge


****************** Get neighbor region for all com ******************
tostring com, replace
replace com = "0"+com if strlen(com) == 4
merge 1:1 com using loc_i_raw.dta 
drop if _merge == 2
drop _merge
destring com, replace


********* Other side ****************
geonear com x_com y_com using loc_j_raw.dta , neighbors(dep x_dep y_dep) nearcount(2)
gen other_dep = nid1 if nid1 != code_dept
replace other_dep = nid2 if nid1 == code_dept

egen home_lab_mkt = sum(population), by(code_dept)
egen other_lab_mkt = sum(population), by(other_dep)

****************** Gen treatment ******************

replace dist_dep = -dist_dep if other_lab_mkt > home_lab_mkt
replace dist_instr = -dist_instr if other_lab_mkt > home_lab_mkt 
gen treat = (other_lab_mkt < home_lab_mkt)

****************** Regressions ******************

keep if abs(dist_dep) <= 25

****************** Dist border ******************
foreach i in dist_river sq3 altitude d_diocese inter_roman cassini_road {
qui: rdrobust `i' dist_dep , c(0) 
rdrobust `i' dist_dep , c(0) fuzzy(dist_inst)
est store ATE_`i'
local N_h = e(N_h_l) + e(N_h_r) 
disp(`N_h')
}

esttab ATE_dist_river ATE_sq3 ATE_altitude ATE_d_diocese ATE_inter_roman ATE_cassini_road, /// 
b(%9.3f) se(%9.3f) stats( N N_h , ///
fmt(3 %9.4gc)) nolines nogaps star(* 0.10 ** 0.05 *** 0.01) ///
keep(RD_Estimate) coef(RD\_Estimate  "ATE") 
